Hybrid NRG-DMRG approach to real-time dynamics of quantum impurity systems 
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A hybrid approach to nonequilibrium dynamics of quantum impurity systems is presented. The 
numerical renormalization group serves as a means to generate a suitable low-energy Hamiltonian, 
allowing for an accurate evaluation of the real-time dynamics of the problem up to exponentially 
long times using primarily the time-adaptive density-matrix renormalization group. We extract the 
decay time of the interaction-enhanced oscillations in the interacting resonant-level model and show 
their quadratic divergence with the interaction strength U. Our numerical analysis is in excellent 
agreement with analytic predictions based on an expansion in 1/U. 
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Introduction. The description of strong electronic 
correlations far from thermal equilibrium poses an enor- 
mous theoretical challenge. At the root of the problem 
lies the nonequilibrium density operator which is not ex- 
plicitly known in the presence of interactions. Of particu- 
lar relevance are quench and driven dynamics realized in 
pump-probe experiments [l|, 0| , atomic traps 0, 0] , and 
nanodevices Q , where the full time evolution of the 
density operator should, in principle, be tracked. 

Quantum impurities systems (QIS) have regained con- 
siderable attention over the past 15 years due to the 
advent of carefully designed nanodevices. These generi- 
cally consist of a few locally interacting degrees of free- 
dom, typically a quantum dot, in contact with macro- 
scopic leads. Since driven dynamics in nanodevices is 
of practical relevance to quantum computing and quan- 
tum control, considerable efforts were mounted in recent 
years toward devising approaches capable of treating the 
nonequilibrium state in QIS. 

Significant analytical progress in the calculation of 
real-time dynamics was achieved using different adap- 
tations of perturbative renormalization-group ideas [7J- 
9]. However, with the exception of Ref. [l(J, these are 
confined to the weak-coupling regime. Numerical meth- 
ods, such as applications of the time-dependent density- 
matrix renormalization group (TD-DMRG) [H E3 to 

QIS OH], 

the time-dependent numerical renormaliza- 
tion group (TD-NRG) [15j, an iterated path- integral ap- 
proach [16[, and different continuous-time Monte Carlo 
simulations (lTlllij . are more flexible in the parameter 



regimes theycan treat, but are either restricted to short 



time scales 




or susceptible to finite-size and 
151 . Indeed, finite-size represen- 



discretization errors [R 
tations are faced with an inherent difficulty of accurately 
representing the continuum limit even on intermediate 
time scales. 

In this paper, we report the extension of a recent hy- 
brid approach [201 ] that overcomes some of the major ob- 



stacles hampering the description of quench dynamics in 
QIS. The basic idea is to exploit the outstanding capa- 
bilities of the TD-NRG to bridge over vastly different 
time scales in order to systematically construct an effec- 
tive low-energy Hamiltonian, whose real-time dynamics 
can be calculated using complementary approaches that 
do not rely on the special structure of the Wilson chain. 
In this manner, one can largely eliminate discretization 
errors inherent to the Wilson chain while boosting the 
complementary approach to times scales orders of mag- 
nitude beyond its natural capabilities. As a proof of prin- 
ciple, we have hybridized in Ref. ^ the TD-NRG with 
the Chebyshev expansion technique (CET) [2l|. In this 
paper, we demonstrate the full power of the approach by 
hybridizing the TD-NRG with the TD-DMRG [l2|. 

Focusing on the interacting resonant-level model 
(IRLM) " 



23[, we show that one can essentially elim- 



inate discretization errors on all time scales of interest 
by constructing a suitable hybrid chain. This, in turn, 
allows for a thorough examination of the interaction- 
enhanced oscillations first reported in Ref. [2(J, yielding 
excellent agreement with analytical predictions for their 
frequency and damping time. The latter is shown to di- 
verge quadratically with the interaction strength, demon- 
strating that relaxation to equilibrium can involve new 
time scales far longer than the thermodynamic ones. 

Hybrid- NRG. We begin with a concise derivation of 
the hybrid-NRG [20|. The Hamiltonian H = Hbath + 
%imp + 'Hmix of a quantum impurity problem consists 
of three parts: Hbath models the continuous bath, T-li m p 
represents the decoupled impurity, and 'Hmix describes 
the coupling between the two subsystems. For t < 0, the 
entire system is assumed to be characterized by a den- 
sity operator po associated with an initial Hamiltonian 
W' . Specifically, po can either be the equilibrium density 
operator corresponding to TV , or may project onto one 
of its low-lying eigenstates, typically the ground state. 
At time t = 0, a static perturbation is abruptly switched 



2 



on such that W —> . Our goal is to track the time 
evolution of local expectation values: 0(t) = Tr{p(t)0} 
with p(t) = e~ imf p e imf . 

In Wilson's numerical rcnormalization group 
(NRG) [24[, Hbath is discretized logarithmically us- 
ing a dimensionless parameter A > 1, and mapped 
onto a semi-infinite chain whose open end is coupled to 
the impurity via H m i x . Wilson's chain is characterized 
by exponentially decreasing hopping matrix elements 
t rn oc Dhr m / 2 , defining a natural separation of scales. 
This enables an iterative diagonalization of H, where at 
each step only the lowest vV s eigenstates are retained. 
Terminating the procedure after N steps, the collection 
of states discarded after each iteration combine to form 
a complete basis set of approximate NRG eigenstates of 
% on the iV-site chain [15|. The expectation value of 
any local operator O can be formally expressed as fl5| 



N trun 



°(t > o) = E E °Zp™,m 



(i) 



where r and s run over the NRG eigenstates of at 
iteration m < N, 0™ s is the matrix representation of O 
at that iteration, and p™ r (t) is the corresponding time- 
dependent reduced density matrix. The restricted sum 
over r and s requires that at least one of these states is 
discarded at iteration m. 

Partitioning the sum over m into m < M and M < m 
at some arbitrary but fixed M < N, Eq. ([T]) is recast as 

M trun 

°(* > °) = E E °r,sP?,M + 01+ p{t)l + M }, 

m— r.s 

( 2 ) 

where 1^ projects onto the subspace retained at the con- 
clusion of iteration M. Equation ([2j is formally exact, 
relying solely on the completeness of our basis set flol ]. 
It has the following interpretation. At each energy scale 
DA"'"/ 2 with m < M, only those terms involving at least 
one discarded high-energy state contribute to 0(t), leav- 
ing the contribution of the low-energy subspace retained 
at the conclusion of iteration M. In the process, the NRG 
has produced an effective quantum-impurity Hamiltonian 
Hm+1 with the reduced bandwidth D e s oc DA~ M / 2 : 



n M+1 =J2 E k'\k;M)(k;M\ 

k 

N-1 

+ E E*™ i ^/m+l,,/™, ! ' + H.C.}i+. (3) 

m=M v 

Here, „ creates an electron of flavor (spin) v on 
the chain site m, \k;M) labels the kept NRG eigen- 
states at iteration M, and E^f are the corresponding 
NRG eigenenergies. Usually, one would proceed with the 
NRG to itcratively diagonalize "Hm+i- Here, we follow 
a different route: (i) we abandon the traditional Wilson 



chain and seek an optimal choice for the hopping ampli- 
tudes t m with m > M; (ii) 7-Lm+i is used as input for 
our complementary method of choice in order to com- 
pute Tr{ll / Ol| / /5(i)l| / }; (iii) employing the standard 
NRG approximation pf r (t) e *(E?-E?)t ^( m ) [lij], 

p™r(M) as produced by our method of choice is feed 
back into the TD-NRG to account for the remaining high- 
energy dynamics in Eq. ^j. 

In this paper, we supplement the TD-NRG with 



the adaptive TD-DMRG [121 ]. The system is assumed 
to initially occupy the ground state |V>o) of H l M+1 , 
constructed using the DMRG @. Accordingly, p(t) 
equals \ip{t))(ilj(t)\ with \ip(t)) = e~ tmf \ip ), result- 
ing in Tr{i+6i+p(t)i+ } = (xM(t)\0\ X M(t)). Here 
|Xm(£)) = IjtflV'W) is the projection of \ip(t)) onto the 
low-energy subspace retained at the conclusion of itera- 
tion M. Although Eq. ^ is formally exact for arbitrary 
M, the larger is N x = (xm(£)|xm (*)) < 1 the smaller 
is the contribution of the sum on the right-hand side of 
Eq. ([2]). If (1 — N x ) < e for some small number e, then 
|Xm(£)) = IV'W) + C(v / e)j an d the major contribution to 
the real-time dynamics originates from {xm ($)\0\xm (fit) ■ 

A proper choice of M is important. Initially, the NRG 
level flows [24j of W and W/ are nearly identical. We 
choose M to be a characteristic iteration after which the 
two level flows begin to significantly deviate from one an- 
other [Ulliil]. The corresponding energy scale, DA~ AI / 2 , 
is typically of order the energy difference between W l and 
W. By that choice of M, the major contribution to 0(t) 
stems from (xM(t)\0\xM(t)) . We approximate |xm(£)) 
by exp(—iHj^ +1 t)i + i \ipo), adopting the NRG philosophy 
that excitations on different energy scales are only weakly 
coupled [20[. Thus, the NRG generates a suitable low- 
energy Hamiltonian, allowing for the real-time dynamics 
of the problem to be explored on the exponentially long 
time scale l/D cS oc k M / 2 using mainly the TD-DMRG. 

The model. We shall demonstrate our hybrid NRG- 
DMRG approach by investigating the quench dynamics 
in the IRLM, defined by the Hamiltonian 



U = y , e k C k c k 



k 
U 



E d dU ■ 
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{c\d + He.} 
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k,k' 



4 c 



k' ■ 



2 



(4) 



Here, d) creates an electron on the impurity level, ct 
creates a band electron with momentum k, N k denotes 
the number of distinct k points, and : c\c k , := c\c k , — 
5k,k'9(— £fc) stands for normal ordering with respect to the 
filled Fermi sea. The basic energy scales in the problem 
include the level energy E4 and the hybridization width 
To = irgV 2 , where g is the conduction-electron density 
of states at the Fermi energy. The effect of the contact 
interaction is to renormalize the hybridization width at 
resonance according to Tq — s- r e g w D(Fq/ 'D) 1 ^ 1+a ' > , 
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FIG. 1. (Color online) Time evolution of rid(t) on a double 
Wilson chain, following a sudden change of the level energy 
from E d to E s d . Here T /D = 1(T 5 and (7 = 0. The full red 
line depicts the exact solution on the double Wilson chain 
(obtained by exact diagonalization), the dashed green line 
displays the hybrid NRG-DMRG, and the full black line is 
the exact analytical continuum-limit solution [l5| in the wide- 
band limit. Chain parameters: M = 29, Ai = 1.8, A2 = 1.02, 
A^nter = 4, and N = 180. N s = 50 states are retained both in 
the NRG and in the course of the TD-DMRG. Inset: Exact 
nd(t) on a pure Wilson chain, for E d — — E d = 10Fo and 
different A's. Here N = 500 (1000) for A > 1.1 (A < 1.05). 



FIG. 2. (Color online) Same as Fig. [H for E s d = -E\ = r e fT 
and different values of U. Fo was adjusted separately for each 
value of U so as to maintain T e g/D = 10~ 5 . All other chain, 
NRG, and TD-DMRG parameters are the same as in Fig. \T\ 



to be the customary Wilson hopping amplitudes [241 ] with 
the discretization parameter Ai. Further down the chain 
a second, smaller discretization parameter A2 is used, 
with a magnitude close to but larger than one [loj]. To 
reduce internal reflections, the transition from Ai to A2 
is smoothed according to ijvf+m = \ 



-1/2 



M+m-l 



with 



with a = 25 - S 2 and 5 = (2/tt) arctan(7T£>C//2) (D is 
the bandwidth). Although the thermodynamics of the 
IRLM was investigated over 30 years ago [22|, H^] , there 
has been a recent surge of interest in its nonequilibrium 
properties, particularly for a biased two-lead setting 0, 
la, iHlil- Focusin g on the single-band version of Eq. Q , 



we consider an abrupt shift of the level energy at time 
with the goal of tracking the time 



t 



from E l d to E J d 



evolution of the level occupancy, Ud(t) = (dJ(t)d(t)). 

Hybrid chain. There are two sources of deviations 
from the continuum limit when considering quench dy- 
namics on a pure Wilson chain: (i) internal reflections of 
currents caused by the exponentiall y d ecreasing hopping 
matrix elements along the chain [20l.l28| (leading, in turn, 
to an exponential slowing down of the transport veloc- 
ity); (ii) reflections at the end of the finite-size chain that 
propagate back to the impurity. While the former source 
of error is eliminated for an ordinary tight-binding chain, 
the latter point is unavoidable in nearly all practical cal- 
culations as the total chain length is limited by computa- 
tional demands. Thus, one would like to simultaneously 
minimize the internal reflections and the transport veloc- 
ity down the chain to accurately access long times. 

Guided by these considerations, we found it advan- 
tageous to use a double Wilson chain, constructed by 
patching two separate Wilson chains. The first M hop- 
ping matrix elements t m with < m < M — 1 arc taken 



Ai-^Mm + 1), 0<m<iV intcr , 



A, 



< TO. 



(5) 



Results. The merit of such a double Wilson chain is 
demonstrated in Fig.[Tl where rid(t) is plotted following a 
quench from E d to E d . The interaction U is set to zero, to 
facilitate comparison with an exact continuum-limit solu- 
tion in the wide-band limit [15| , as well as with an exact 
numerical solution on the hybrid chain, obtained by ex- 
act diagonalization of the single-particle eigenmodes. We 
set To/D = 10~ 5 , placing the basic time scale 1/To or- 
ders of magnitude beyond the reach of pure TD-DMRG. 
The parameter M was chosen such that DA 1 M ^ 2 w 20To 
is twice the maximal value of \Ed\ used. 

Evidently, deviations between the continuum limit, the 
exact solution on the hybrid chain, and the hybrid NRG- 
DMRG approach are hardly discernible up to long time 
scales, well after the occupancy has relaxed to its new 



equilibrium value. For E t 
agreement persists up to t 
three curves begin to separate. For E d 



> 



30/r 



To, the excellent 
at which point all 
-E\ - 10T , 



the agreement extends up to slightly above 15/ro. As 
analyzed in the inset, an impractically small discretiza- 
tion parameter A rj A 2 = 1.02 is needed to achieve a 
comparable representation of the continuum limit using 
a pure Wilson chain. 

In Fig. [21 we analyze the effect of a finite U on rid(t). 
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FIG. 3. (Color online) (a) Interaction-enhanced oscillations 
(red line) for E f d = -E\ = 3r cff and U/D = 10, along 
with a fit to Eq. {7} (black crosses), (b) Same as (a), for 

e ! , = -El = r c 



E f d /r cB = -Ei/r eB 



(c) The fitted frequency f2 vs U for 
= 1 (black circles), 3 (red circles), and 
10 (green circles). Full lines display the analytical strong- 
coupling expression for Q, employing the numerical values 
of V. Arrows on the right-hand side mark the asymptotic 
U — > oo values of Q. (d) Same as (c) for the decay time r. 
All remaining parameters are the same as in Fig. [2] 



At low energies, the IRLM is equivalent to its noninter- 
acting counterpart, both describing a phase-shifted Fermi 
liquid. Near resonance, the effect of U in equilibrium is 
to renormalizc Tq to T e g, hence one may expect nd(t) 
to follow the same curves as in Fig. [T] upon substituting 
To — > r c ff. This, however, is not the case. In Fig. [5] 
we adjusted Tq separately for each value of U so as to 
maintain a fixed T e s/D = 10 -5 . To this end, we fixed 
Ed/D = — 10~ 5 and scanned r using the hybrid NRG- 
DMRG until a ground-state occupancy of n c i = 0.75 
was reached. As first reported in Ref. [201 ] . interaction- 
enhanced oscillations gradually develop in nd(t) upon 
increasing U. For large U, exemplified by U/D = 5, 
these oscillations decay on a time scale much longer than 
l/r e ff, revealing the emergence of a new time scale un- 
related to the thermodynamic ones. Note that, as for 
U = 0, our curves appear to faithfully represent the con- 
tinuum limit up to t > 30/r c ff for this moderate quench. 

The extended time scales and supreme accuracy of 
the hybrid NRG-DMRG allow for a detailed quan- 
titative analysis of the interaction-enhanced oscilla- 
tions, which was previously impossible using the hybrid 
NRG-Chebyshev approach [20|. The understanding of 
the interaction-enhanced oscillations employs a strong- 
coupling expansion in 1/U. For U — > oo, the impu- 
rity level and zeroth Wilson shell decouple from the rest 
of the chain, being confined to a combined valence of 



one. Within this subspace, the two eigenstates of the 
local Hamiltonian have the energies e± = (Ed/ 2) ± 
\J (E d /2) 2 + V 2 , hence n d {t) displays quantum beats 
with the frequency Q = e+ - e_ = 2yJ(E d /2) 2 + V 2 
For large but finite U, the coherent oscillations are 
damped by the residual coupling to the rest of the chain, 
which introduces a finite lifetime of the state e+. We 
expand to order 1/U about the U — > oo limit, which 
yields the residual coupling to the rest of the chain. Us- 
ing Fermi's golden rule, the decay time is found to be 



8D 



V 2 



* 2 U) ^{E d /2) 2 + V 2 



(6) 



Further neglecting rearrangements of the bath electrons 
(themselves being controlled by 1/U), we deduce the 
functional form 



e cos (fit) VI 



-t/2r, 



n d (t) = A 

+ n cq (l - e-*/ T ) +n d (0)e 



e~*/ T cos 2 i 

-t/r 



e~ t/T smt 



(7) 



where n cq is the equilibrated long-time occupancy, while 
A and 9 have no direct relation to any simple observable. 
Panels (a) and (b) of Fig.|3]show typical fits of rid(t) to 



Eq. ©, using r, fi, A, 6, and n rf (0) = 1 



as fitting- 



parameters. The fitting range, t ■ T c g < 10, was carefully 
chosen to exclude any discretization error. Evidently, 
Eq. |7]) well describes the numerical curves, further vali- 
dating the expansion in 1/U. The extracted values of fl 
and r, plotted in panels (c) and (d) for different quenches, 
practically coincide above U/D w 10 with the analytical 
predictions, confirming, in particular, that t (xU 2 . 

Summary. A new hybrid NRG-DMRG approach was 
devised that largely eliminates discretization errors ham- 
pering the TD-NRG, while boosting the TD-DMRG to 
time scales orders of magnitude beyond its natural reach. 
The approach allows access to exceptionally long times 
with unparallcl accuracy, as demonstrated by a detailed 
analysis of the interaction-enhanced oscillations in the 
IRLM. These outstanding capabilities open the door, so 
we hope, to accurate investigations of systems and cou- 
pling regimes that so far remained well beyond reach. 

Acknowledgments This work was supported by the 
German-Israeli Foundation through grant no. 1035-36.14, 
and by the Deutsche Forschungsgemeinschaft under AN 
275/6-2 (F.G. and F.B.A). 



[1] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, 
H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, and 
M. Wolf, Phys. Rev. Lett. 97, 067402 (2006). 

[2] F. Schmitt, P. S. Kirchmann, U. Bovensiepen, R. G. 
Moore, L. Rettig, M. Krenz, J.-H. Chu, N. Ru, L. Per- 
fetti, D. H. Lu, M. Wolf, I. R. Fisher, and Z.-X. Shen, 
Science 321, 1649 (2008). 



5 



[3] M. Greiner, O. Mandel, T. W. Hansch, and I. Bloch, 

Nature 419, 51 (2002); T. Kinoshita, T. Wenger, and D. 

S. Weiss, Nature 440, 900 (2006). 
[4] S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. 

Schollwock, J. Eisert, and I. Bloch, Nature Physics 8, 

325 (2012). 

[5] J. M. Elzerman, R. Hanson, L. H. W. van Beveren, B. 
Witkamp, L. M. K. Vandersypen, and L. P. Kouwen- 
hoven, Nature 430, 431 (2004). 
[6] J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. 
Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and 
A. C. Gossard, Science 309, 2180 (2005). 
[7] M. Pletyukhov, D. Schuricht, and H. Schoeller, Phys. 
Rev. Lett. 104, 106801 (2010); S. Andergassen, M. Ple- 
tyukhov, D. Schuricht, H. Schoeller, and L. Borda, Phys. 
Rev. B 83, 205103 (2011); C. Karrasch, S. Andergassen, 
M. Pletyukhov, D. Schuricht, L. Borda, V. Meden, and 
H. Schoeller, Europhys. Lett. 90, 30003 (2010). 
[8] A. Hackl and S. Kehrein, J. Phys. C 21, 015601 (2009); 
A. Hackl, M. Vojta, and S. Kehrein, Phys. Rev. B 80, 
195117 (2009); P. Wang and S. Kehrein, Phys. Rev. B 
82, 125124 (2010). 
[9] D. M. Ken nes, S. G. Jakobs, C. Karrasch, and V. Meden, 
Report no. larXiv:1111.6982l 

[10] M. Pletyukhov and H. Schoeller, Report no. 
larXiv:1201. 62951 

[11] G. Vidal, Phys. Rev. Lett. 93, 040502 (2004); A. J. Da- 
ley, C. Kollath, U. Schollwock, and G. Vidal, J. Stat. 
Mech.: Theor. Exp. P04005 (2004); S. R. White and A. 
E. Fciguin, Phys. Rev. Lett. 93, 076401 (2004). 

[12] For a recent review on the DMRG and TD-DMRG, see 
U. Schollwock, Ann. Phys. 326, 96 (2011). 

[13] P. Schmitteckert, Phys. Rev. B 70, 121302 (2004); A. 
Branschaedel, G. Schneider, and P. Schmitteckert, Ann. 
Phys. 522, 657 (2010). 

[14] K. A. Al-Hassanieh, A. E. Feiguin, J. A. Riera, C. A. 
Biisser, and E. Dagotto, Phys. Rev. B 73, 195304 (2006); 
L. G. G. V. Dias da Silva, F. Heidrich-Meisner, A. E. 
Feiguin, C. A. Biisser, G. B. Martins, E. V. Anda, and 
E. Dagotto, Phys. Rev. B 78, 195317 (2008); F. Heidrich- 



Meisner, A. E. Feiguin, and E. Dagotto, Phys. Rev. B 79, 
235336 (2009). 

[15] F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 

(2005) ; Phys. Rev. B 74, 245113 (2006). 

[16] S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys. 

Rev. B 77, 195316 (2008). 
[17] L. Miihlbacher and E. Rabani, Phys. Rev. Lett. 100, 

176403 (2008). 

[18] P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B 79, 
035320 (2009); P. Werner, T. Oka, M. Eckstein, and A. 
J. Millis, Phys. Rev. B 81, 035108 (2010). 

[19] M. Schiro and M. Fabrizio, Phys. Rev. B 79, 153302 
(2009). 

[20] E. Eidelstein, A. Schiller, F. Giittge, and F. B. Anders, 

Phys. Rev. B 85, 075118 (2012). 
[21] H. Tal Ezer and R. Kosloff, J. Chem. Phys. 81, 3967 

(1984). 

[22] P. W. Vigman and A. M. Finkelstein, Zh. Eksp. Theor. 

Fiz. 75, 204 (1978) [Sov. Phys. JETP 48, 102 (1978)]. 
[23] P. Schlottmann, Phys. Rev. B 22, 622 (1980). 
[24] K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975); R. Bulla, 

T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 

(2008) . 

[25] S. R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. 

Rev. B 48, 10345 (1993). 
[26] P. Mehta and N. Andrei, Phys. Rev. Lett. 96, 216802 

(2006) . 

[27] E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. 

Lett. 101, 140601 (2008). 
[28] P. Schmitteckert, J. Phys.: Conference Series 220, 

012022 (2010). 

[29] A. Weichselbaum, F. Verstraete, U. Schollwock, J. I. 
Cirac, and Jan von Delft, Phys. Rev. B 80, 165117 

(2009) . 

[30] Another motivation for u sing a modified Wilson chain 
was recently given in Ref. [23]. These authors adopted a 
different point of view and a different construction based 
upon Wilson's original discretization of the bath. 



